Morphology.f90 Source File
Module to deal with river and basin morphology
Source Code
!! Module to deal with river and basin morphology !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.2 - 22nd April 2024 ! ! | version | date | comment | ! |----------|--------------|----------| ! | 0.1 | 17/Aug/2009 | Original code | ! | 0.2 | 23/Oct/2010 | Compute flow accumulation grid | ! | 0.3 | 04/Dec/2012 | Compute longest flow length | ! | 0.4 | 03/Apr/2013 | Derive aspect subroutine | ! | 0.5 | 10/Apr/2013 | Derive curvature subroutine | ! | 0.6 | 13/Apr/2013 | Fixed DeriveSlope to match ArcGis algorithm | ! | 0.7 | 16/Apr/2013 | DrainageDensity rooutine | ! | 0.8 | 03/May/2013 | TopographicIndex routine | ! | 0.9 | 10/Jun/2013 | CurvatureMicroMet routine | ! | 1.0 | 09/Oct/2013 | LongestFlowPathSlope | ! | 1.1 | 03/Feb/2021 | Compute Centroid | ! | 1.2 | 22/Apr/2024 | ESRI and GRASS flow direction can be used | ! ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! ! This file is part of ! ! MOSAICO -- MOdular library for raSter bAsed hydrologIcal appliCatiOn. ! ! Copyright (C) 2011 Giovanni Ravazzani ! !### Code Description ! Language: Fortran 90. ! ! Software Standards: "European Standards for Writing and ! Documenting Exchangeable Fortran 90 Code". ! !### Module Description ! library to deal with river and basin morphology MODULE Morphology ! Modules used: USE DataTypeSizes, ONLY: & !Imported parameters: short, & long, & float USE GridLib, ONLY: & !Imported type definitions: grid_integer, & grid_real, & ! Imported routines: NewGrid, & GridDestroy, & ExportGrid, & !Imported parameters: ESRI_ASCII, & ESRI_BINARY, & NET_CDF USE GridOperations, ONLY: & !Imported routines: IsOutOfGrid, & GetIJ, & CellArea, & GetXY, & GridResample, & ExtractBorder, & CellLength, & ASSIGNMENT(=) USE GeoLib, ONLY: & !Imported definitions: Coordinate, & !Imported variables: point1, & point2, & !Imported routines: Distance, & ASSIGNMENT (=) USE LogLib, ONLY: & !Imported routines: Catch USE Units, ONLY: & !Imported parameters: pi USE StringManipulation, ONLY: & !Imported routines: ToString, & StringToLower IMPLICIT NONE PUBLIC :: HortonOrders PUBLIC :: CheckOutlet PUBLIC :: DownstreamCell PUBLIC :: FlowAccumulation PUBLIC :: DeriveSlope PUBLIC :: DeriveAspect PUBLIC :: DeriveCurvature PUBLIC :: CurvatureMicroMet PUBLIC :: BasinDelineate PUBLIC :: CellIsSpring PUBLIC :: LongestFlowLength PUBLIC :: LongestFlowPathSlope PUBLIC :: Perimeter PUBLIC :: Centroid PUBLIC :: SetFlowDirectionConvention PRIVATE :: ConfluenceIsAround PRIVATE :: BasinMask PRIVATE :: BasinArea PRIVATE :: CentroidGridReal PRIVATE :: CentroidGridInteger INTERFACE Centroid MODULE PROCEDURE CentroidGridReal MODULE PROCEDURE CentroidGridInteger END INTERFACE !flow direction INTEGER (KIND = short) :: NW !! North-West INTEGER (KIND = short) :: N !! North INTEGER (KIND = short) :: NE !! North-East INTEGER (KIND = short) :: W !! West INTEGER (KIND = short) :: E !! East INTEGER (KIND = short) :: SW !! South-West INTEGER (KIND = short) :: S !! South INTEGER (KIND = short) :: SE !! South-East !flow direction conventions INTEGER (KIND = short), PARAMETER :: ESRI = 1, GRASS = 2 INTEGER (KIND = short) :: flowDirectionConvention !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Set flow direction convention: ! ! - ESRI: NW=32, N=64, NE=128, W=16, E=1, SW=8, S=4, SE=2 ! - GRASS: NW=3, N=2, NE=1, W=4, E=8, SW=5, S=6, SE=7 ! SUBROUTINE SetFlowDirectionConvention & ! (string) IMPLICIT NONE !Arguments with intent in: CHARACTER (LEN = *), INTENT (IN) :: string !-------------------------------end of declarations---------------------------- SELECT CASE ( TRIM (StringToLower(string) ) ) CASE ('esri') flowDirectionConvention = ESRI NW = 32; N = 64; NE = 128; W = 16; E = 1; SW = 8; S = 4; SE = 2 CASE ('grass') flowDirectionConvention = GRASS NW = 3; N = 2; NE = 1; W = 4; E = 8; SW = 5; S = 6; SE = 7 CASE DEFAULT CALL Catch ('error', 'Morphology', 'flow direction convention not found: ', & argument = TRIM (string) ) END SELECT RETURN END SUBROUTINE SetFlowDirectionConvention !============================================================================== !| Description: ! returns the coordinate of the centroid given a `grid_integer` mask, ! in the same coordinate reference system of `grid_integer` SUBROUTINE CentroidGridInteger & ! (grid, ccoord) IMPLICIT NONE !Arguments with intent in: TYPE (grid_integer), INTENT (IN) :: grid !!input grid !Arguments with intent out or inout TYPE (Coordinate), INTENT (OUT) :: ccoord !! coordinate of computed centroid !local declarations: INTEGER (KIND = short) :: cj, ci !!coordinate INTEGER (KIND = short) :: i, j INTEGER (KIND = short) :: countcell !--------------------------------end of declaration---------------------------- !set coordinate reference system of centroid ccoord % system = grid % grid_mapping cj = 0 ci = 0 countcell = 0 DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN countcell = countcell + 1 cj = cj + j ci = ci + i END IF END DO END DO IF ( countcell > 0) THEN cj = cj / countcell ci = ci / countcell ELSE CALL Catch ('error', 'Morphology', 'area zero found while computing centroid coordinates') END IF CALL GetXY (ci, cj, grid, ccoord % easting, ccoord % northing ) RETURN END SUBROUTINE CentroidGridInteger !============================================================================== !| Description: ! returns the coordinate of the centroid given a `grid_integer` mask, ! in the same coordinate reference system of `grid_integer` SUBROUTINE CentroidGridReal & ! (grid, ccoord) IMPLICIT NONE !Arguments with intent in: TYPE (grid_real), INTENT (IN) :: grid !!input grid !Arguments with intent out or inout TYPE (Coordinate), INTENT (OUT) :: ccoord !! coordinate of computed centroid !local declarations: INTEGER (KIND = short) :: cj, ci !!coordinate INTEGER (KIND = short) :: i, j INTEGER (KIND = short) :: countcell !--------------------------------end of declaration---------------------------- !set coordinate reference system to centroid ccoord % system = grid % grid_mapping cj = 0 ci = 0 countcell = 0 DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN countcell = countcell + 1 cj = cj + j ci = ci + i END IF END DO END DO IF ( countcell > 0) THEN cj = cj / countcell ci = ci / countcell ELSE CALL Catch ('error', 'Morphology', 'area zero found while computing centroid coordinates') END IF CALL GetXY (ci, cj, grid, ccoord % easting, ccoord % northing ) RETURN END SUBROUTINE CentroidGridReal !============================================================================== !| Description: ! returns a grid_integer containing Horton orders. Horton orders are ! computed on the entire space-filled basin. SUBROUTINE HortonOrders & ! (flowDirection,orders,basinOrder) IMPLICIT NONE !Arguments with intent in: TYPE (grid_integer), INTENT (IN) :: flowDirection !Arguments with intent out or inout TYPE (grid_integer), INTENT (INOUT) :: orders INTEGER, OPTIONAL, INTENT (OUT) :: basinOrder !!the maximum order of the basin !local declarations: LOGICAL :: confluence !!true if confluence LOGICAL :: outlet !!true if basin outlet INTEGER :: row, col !!current cell INTEGER :: iDown, jDown !!downstream cell INTEGER :: numConf !!number of confluences INTEGER :: order !! Horton order INTEGER :: cellsCount INTEGER :: i, j !--------------------------------end of declaration---------------------------- order = 1 numConf = 1 DO WHILE (numConf > 0) ! se non trovo confluenze ! di classe order ! l'operazione è terminata CALL Catch ('info', 'Morphology', 'Elaborating reaches of stream order: ', & argument = ToString(order)) numConf = 0 !-----follow the reach till a confluence or a basin outlet------ DO j = 1,orders % jdim DO i = 1,orders % idim IF(CellIsSpring(i,j,flowDirection)) THEN !found a spring row = i col = j outlet = .FALSE. confluence = .FALSE. cellsCount = 0 orders % mat(i,j) = 1 DO WHILE (.NOT. outlet) ! follow the reach till the basin outlet IF (orders % mat(row,col) == order ) THEN cellsCount = cellsCount + 1 ENDIF CALL DownstreamCell(row, col, & flowDirection%mat(row,col), & iDown, jDown) IF (cellsCount >= 1 ) THEN !I am in the reach of that order !check if downstream cell is a confluence to increment horton order !Downstream the confluence, till the basin outlet, as temptative value, !order is increased by 1 (order + 1) IF ( .NOT. confluence ) THEN CALL ConfluenceIsAround (iDown, jDown, row, col, & flowDirection,confluence,orders,order) IF(confluence) numConf = numConf + 1 ENDIF outlet = CheckOutlet (row,col,iDown,jDown,flowDirection) IF (.NOT. outlet) THEN IF (.NOT. confluence) THEN orders % mat(iDown,jDown) = order ELSE orders % mat(iDown,jDown) = order + 1 ENDIF ENDIF ENDIF ! cellsCount >= 1 outlet = CheckOutlet(row,col,iDown,jDown,flowDirection) !loop row = iDown col = jDown END DO ENDIF ENDDO ENDDO !ciclo sulla matrice ordini !------------------------------------------------------------------------------ order = order + 1 ENDDO IF ( PRESENT (basinOrder) ) THEN basinOrder = order - 1 END IF END SUBROUTINE HortonOrders !============================================================================== !| Description: ! Scan the eigth cells surrounding the center cell (row,col) (neglecting ! the upstream cell (i,j) ) to find if a confluence is present. ! The confluence cell must be of the same order or not yet defined. SUBROUTINE ConfluenceIsAround & ! (row,col,i,j,flowDir,confluence,orders,Norder) IMPLICIT NONE TYPE (grid_integer), INTENT(IN):: flowDir TYPE (grid_integer), INTENT(IN):: orders INTEGER, INTENT(IN):: i,j,row,col,Norder LOGICAL, INTENT(OUT):: confluence !----------------------------------------------------------- IF(.NOT. IsOutOfGrid(row,col+1,flowDir) ) THEN IF(flowDir%mat(row,col+1) == W .AND.& (row /= i .OR. col+1 /= j).AND. & (orders%mat(row,col+1) == orders%nodata.OR. & orders%mat(row,col+1) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row+1,col+1,flowDir) ) THEN IF(flowDir%mat(row+1,col+1) == NW .AND.& (row+1 /=i .OR. col+1 /= j).AND. & (orders%mat(row+1,col+1) == orders%nodata.OR. & orders%mat(row+1,col+1) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row+1,col,flowDir) ) THEN IF(flowDir%mat(row+1,col) == N .AND.& (row+1 /=i .OR. col /= j).AND. & (orders%mat(row+1,col) == orders%nodata.OR. & orders%mat(row+1,col) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row+1,col-1,flowDir) ) THEN IF(flowDir%mat(row+1,col-1) == NE .AND.& (row+1 /= i .OR. col-1 /= j).AND. & (orders%mat(row+1,col-1) == orders%nodata.OR. & orders%mat(row+1,col-1) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row,col-1,flowDir) ) THEN IF(flowDir%mat(row,col-1) == E .AND.& (row /= i .OR. col-1 /= j).AND. & (orders%mat(row,col-1) == orders%nodata.OR. & orders%mat(row,col-1) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row-1,col-1,flowDir) ) THEN IF(flowDir%mat(row-1,col-1) == SE .AND.& (row-1 /= i .OR. col-1 /= j).AND. & (orders%mat(row-1,col-1) == orders%nodata.OR. & orders%mat(row-1,col-1) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row-1,col,flowDir) ) THEN IF(flowDir%mat(row-1,col) == S .AND.& (row-1 /= i .OR. col /= j).AND. & (orders%mat(row-1,col) == orders%nodata.OR. & orders%mat(row-1,col) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF IF(.NOT. IsOutOfGrid(row-1,col+1,flowDir) ) THEN IF(flowDir%mat(row-1,col+1) == SW .AND.& (row-1 /= i .OR. col+1 /= j).AND. & (orders%mat(row-1,col+1) == orders%nodata.OR. & orders%mat(row-1,col+1) == Norder) ) THEN confluence = .TRUE. ENDIF ENDIF RETURN END SUBROUTINE ConfluenceIsAround !============================================================================== !| Description: ! if the downstream cell is a nodata value or out of grid space, ! or points toward the current cell, the current cell ! is the basin outlet. LOGICAL FUNCTION CheckOutlet & ! (icurrent, jcurrent, iDown, jDown, flowDirection) IMPLICIT NONE INTEGER, INTENT(in) :: icurrent INTEGER, INTENT(in) :: jcurrent INTEGER, INTENT(in) :: iDown !riga cella di valle INTEGER, INTENT(in) :: jDown !colonna cella di valle TYPE(grid_integer), INTENT(in):: flowDirection !local declarations: INTEGER :: iback, jback !------------------------------end of declaration ----------------------------- CheckOutlet = .FALSE. !first case IF ( IsOutOfGrid (iDown,jDown,flowDirection) ) THEN CheckOutlet = .TRUE. RETURN END IF !second case IF ( flowDirection % mat(iDown,jDown) == flowDirection % nodata) THEN CheckOutlet = .TRUE. RETURN END IF !third case CALL DownstreamCell (iDown,jDown,flowDirection%mat(iDown,jDown),iback,jback) IF ( iback == icurrent .AND. jback == jcurrent) THEN CheckOutlet = .TRUE. RETURN END IF RETURN END FUNCTION CheckOutlet !============================================================================== !| Description: ! returns the position (is,js) of the downstream cell and, optionally, ! the flow path length, considering cardinal and diagonal direction SUBROUTINE DownstreamCell & ! (iin, jin, dir, is, js, dx, grid) IMPLICIT NONE !Arguments with intent in: INTEGER (KIND = short), INTENT (IN) :: iin, jin !!current cell INTEGER (KIND = long), INTENT (IN) :: dir !!flow direction TYPE(grid_integer), INTENT (IN),OPTIONAL:: grid !!used to define coordinate reference system !Arguments with intent out: INTEGER (KIND = short), INTENT (OUT) :: is,js !!downstream cell REAL (KIND = float), INTENT (OUT) ,OPTIONAL :: dx !!flow path length [m] !Local declarations REAL (KIND = float) :: ddx,x,y,xs,ys !--------------------end of declarations--------------------------------------- !east IF ( dir == E ) THEN js = jin + 1 is = iin END IF !south east IF ( dir == SE ) THEN js = jin + 1 is = iin + 1 END IF !south IF ( dir == S ) THEN js = jin is = iin + 1 END IF !south west IF ( dir == SW ) THEN js = jin - 1 is = iin + 1 END IF !west IF ( dir == W ) THEN js = jin - 1 is = iin END IF !north west IF ( dir == NW ) THEN js = jin - 1 is = iin - 1 END IF !north IF ( dir == N ) THEN js = jin is = iin - 1 END IF !north easth IF ( dir == NE ) THEN js = jin + 1 is = iin - 1 END IF IF ((PRESENT(dx)).and.(PRESENT(grid))) THEN CALL GetXY (iin,jin,grid,x,y) CALL GetXY (is,js,grid,xs,ys) point1 % system = grid % grid_mapping point2 % system = grid % grid_mapping point1 % northing = y point1 % easting = x point2 % northing = ys point2 % easting = xs ddx = distance(point1,point2) dx = ddx RETURN ELSE IF ((PRESENT(dx)) .AND. .NOT. (PRESENT(grid))) THEN CALL Catch ('error', 'Morphology', 'missing grid while calculating downstream distance') END IF END SUBROUTINE DownstreamCell !============================================================================== !| Description: ! returns the position (is,js) of the upstream cell and, optionally, ! the flow path length, considering cardinal and diagonal direction SUBROUTINE UpstreamCell & ! (i, j, flowdir, flowacc, is, js, found, dx) IMPLICIT NONE !Arguments with intent in: INTEGER (KIND = short), INTENT (IN) :: i, j !!current cell TYPE (grid_integer), INTENT (IN) :: flowdir !!flow direction map TYPE (grid_real), INTENT (IN) :: flowacc !! flow accumulation map !Arguments with intent out: INTEGER (KIND = short), INTENT (OUT) :: is,js !!upstream cell LOGICAL, INTENT(OUT) :: found !! upsteam cell found REAL (KIND = float), INTENT (OUT) ,OPTIONAL :: dx !!flow path length [m] !Local declarations REAL (KIND = float) :: x,y,xs,ys INTEGER (KIND = short) :: row, col REAL (KIND = float) :: minDeltaArea, deltaArea !--------------------end of declarations--------------------------------------- minDeltaArea = 1000000000. is = 0 js = 0 found = .FALSE. !east cell row = i col = j + 1 IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == W ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !south-east cell row = i + 1 col = j + 1 IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == NW ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !south cell row = i + 1 col = j IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == N ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !south-west cell row = i + 1 col = j - 1 IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == NE ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !west cell row = i col = j - 1 IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == E ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !north-west cell row = i - 1 col = j - 1 IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == SE ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !north cell row = i - 1 col = j IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == S ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF !north-east cell row = i - 1 col = j + 1 IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN IF(flowDir%mat(row,col) == SW ) THEN !check delta area deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col) IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN found = .TRUE. minDeltaArea = deltaArea is = row js = col END IF ENDIF ENDIF IF (found) THEN IF (PRESENT (dx) ) THEN CALL GetXY (i,j,flowDir,x,y) CALL GetXY (is,js,flowDir,xs,ys) point1 % system = flowDir % grid_mapping point2 % system = flowDir % grid_mapping point1 % northing = y point1 % easting = x point2 % northing = ys point2 % easting = xs dx = distance(point1,point2) END IF END IF END SUBROUTINE UpstreamCell !============================================================================== !| Description: ! For each cell, the routine calculates the maximum rate of change in value ! from that cell to its neighbors. Basically, the maximum change in ! elevation over the distance between the cell and its eight neighbors ! identifies the steepest downhill descent from the cell. ! Conceptually, the routine fits a plane to the z-values of a 3 x 3 cell ! neighborhood around the processing or center cell. The slope value of this ! plane is calculated using the average maximum technique. The direction ! the plane faces is the aspect for the processing cell. ! If there is a cell location in the neighborhood with a NoData z-value, ! the z-value of the center cell will be assigned to the location. ! ! References: ! ! Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical ! Information Systems (Oxford University Press, New York), 190 pp. SUBROUTINE DeriveSlope & ! (dem, slope) IMPLICIT NONE !arguments with intent in TYPE(grid_real), INTENT(in):: dem !arguments with intent out TYPE (grid_real), INTENT (out) :: slope ![rad] !local variables INTEGER :: ii,jj REAL (KIND = float) :: dZdX !!rate of change in the x direction REAL (KIND = float) :: dZdY !!rate of change in the y direction REAL (KIND = float) :: a !!elevation of N-W cell REAL (KIND = float) :: b !!elevation of N cell REAL (KIND = float) :: c !!elevation of N-E cell REAL (KIND = float) :: d !!elevation of W cell REAL (KIND = float) :: f !!elevation of E cell REAL (KIND = float) :: g !!elevation of S-W cell REAL (KIND = float) :: h !!elevation of S cell REAL (KIND = float) :: i !!elevation of S-E cell REAL (KIND = float) :: L !!length !------------------------------end of declaration ----------------------------- !allocate slope grid CALL NewGrid (slope,dem) !loop through dem grid idim: DO ii = 1, dem % idim jdim: DO jj = 1, dem % jdim IF (dem % mat(ii,jj) /= dem % nodata) THEN !set cell size L = CellArea (dem,ii,jj) ** 0.5 !north-west cell IF ( .NOT. IsOutOfGrid (ii-1,jj-1,dem) ) THEN IF ( dem % mat (ii-1,jj-1) /= dem % nodata) THEN a = dem % mat (ii-1,jj-1) ELSE a = dem % mat (ii,jj) END IF ELSE a = dem % mat (ii,jj) END IF !north cell IF ( .NOT. IsOutOfGrid (ii-1,jj,dem) ) THEN IF ( dem % mat (ii-1,jj) /= dem % nodata) THEN b = dem % mat (ii-1,jj) ELSE b = dem % mat (ii,jj) END IF ELSE b = dem % mat (ii,jj) END IF !north-east cell IF ( .NOT. IsOutOfGrid (ii-1,jj+1,dem) ) THEN IF ( dem % mat (ii-1,jj+1) /= dem % nodata) THEN c = dem % mat (ii-1,jj+1) ELSE c = dem % mat (ii,jj) END IF ELSE c = dem % mat (ii,jj) END IF !west cell IF ( .NOT. IsOutOfGrid (ii,jj-1,dem) ) THEN IF ( dem % mat (ii,jj-1) /= dem % nodata) THEN d = dem % mat (ii,jj-1) ELSE d = dem % mat (ii,jj) END IF ELSE d = dem % mat (ii,jj) END IF !east cell IF ( .NOT. IsOutOfGrid (ii,jj+1,dem) ) THEN IF ( dem % mat (ii,jj+1) /= dem % nodata) THEN f = dem % mat (ii,jj+1) ELSE f = dem % mat (ii,jj) END IF ELSE f = dem % mat (ii,jj) END IF !south-west cell IF ( .NOT. IsOutOfGrid (ii+1,jj-1,dem) ) THEN IF ( dem % mat (ii+1,jj-1) /= dem % nodata) THEN g = dem % mat (ii+1,jj-1) ELSE g = dem % mat (ii,jj) END IF ELSE g = dem % mat (ii,jj) END IF !south cell IF ( .NOT. IsOutOfGrid (ii+1,jj,dem) ) THEN IF ( dem % mat (ii+1,jj) /= dem % nodata) THEN h = dem % mat (ii+1,jj) ELSE h = dem % mat (ii,jj) END IF ELSE h = dem % mat (ii,jj) END IF !south-east cell IF ( .NOT. IsOutOfGrid (ii+1,jj+1,dem) ) THEN IF ( dem % mat (ii+1,jj+1) /= dem % nodata) THEN i = dem % mat (ii+1,jj+1) ELSE i = dem % mat (ii,jj) END IF ELSE i = dem % mat (ii,jj) END IF !compute the rate of change in the x direction dZdX = ((c + 2.*f + i) - (a + 2.*d + g)) / (8. * L) !compute the rate of change in the y direction dZdY = ((g + 2.*h + i) - (a + 2.*b + c)) / (8. * L) !compute slope slope % mat (ii,jj) = ATAN ( (dZdX ** 2. + dZdY ** 2.) ** 0.5 ) END IF END DO jdim END DO idim RETURN END SUBROUTINE DeriveSlope !============================================================================== !| Description: ! Aspect identifies the downslope direction of the maximum rate of change ! in value from each cell to its neighbors. ! Aspect can be thought of as the slope direction. The values of the output ! raster will be the compass direction of the aspect. ! A moving window visits each cell in the input raster and for each cell ! in the center of the window, an aspect value is calculated using an ! algorithm that incorporates the values of the cell's eight neighbors. ! If any neighborhood cells are NoData, they are first assigned the value ! of the center cell, then the aspect is computed. ! The cells are identified as letters 'a' to 'i', with 'e' representing the ! cell for which the aspect is being calculated. ! The rate of change in the x direction for cell 'e' is calculated with ! the following algorithm: !``` ! [dz/dx] = ((c + 2f + i) - (a + 2d + g)) / 8 !``` ! The rate of change in the y direction for cell 'e' is calculated ! with the following algorithm: !``` ! [dz/dy] = ((g + 2h + i) - (a + 2b + c)) / 8 !``` ! Taking the rate of change in both the x and y direction for ! cell 'e', aspect is calculated using: !``` ! aspect = 57.29578 * atan2 ([dz/dy], -[dz/dx]) !``` ! The aspect value is then converted to compass direction values ! (0–2pi rad), according to the following rule: !``` ! if aspect < 0 ! ! cell = pi/2 - aspect ! ! else if aspect > pi/2 ! ! cell = 2pi - aspect + pi/2 ! ! else ! ! cell = pi/2 - aspect !``` ! Flat areas in the input raster—areas where the slope is zero are assigned ! an aspect of -1. ! ! Reference: ! ! Burrough, P. A. and McDonell, R.A., 1998. Principles of Geographical ! Information Systems (Oxford University Press, New York), p. 190 SUBROUTINE DeriveAspect & ! (dem, aspect) IMPLICIT NONE !arguments with intent in TYPE(grid_real), INTENT(in):: dem !arguments with intent out TYPE (grid_real), INTENT (out) :: aspect !(rad) !local variables INTEGER :: ii,jj REAL (KIND = float) :: dZdX !!rate of change in the x direction REAL (KIND = float) :: dZdY !!rate of change in the y direction REAL (KIND = float) :: aspectIJ REAL (KIND = float) :: a !!elevation of N-W cell REAL (KIND = float) :: b !!elevation of N cell REAL (KIND = float) :: c !!elevation of N-E cell REAL (KIND = float) :: d !!elevation of W cell REAL (KIND = float) :: f !!elevation of E cell REAL (KIND = float) :: g !!elevation of S-W cell REAL (KIND = float) :: h !!elevation of S cell REAL (KIND = float) :: i !!elevation of S-E cell !------------------------------end of declaration ----------------------------- !allocate aspect grid CALL NewGrid (aspect,dem) !loop through dem grid idim: DO ii = 1, dem % idim jdim: DO jj = 1, dem % jdim IF (dem % mat(ii,jj) /= dem % nodata) THEN !north-west cell IF ( .NOT. IsOutOfGrid (ii-1,jj-1,dem) ) THEN IF ( dem % mat (ii-1,jj-1) /= dem % nodata) THEN a = dem % mat (ii-1,jj-1) ELSE a = dem % mat (ii,jj) END IF ELSE a = dem % mat (ii,jj) END IF !north cell IF ( .NOT. IsOutOfGrid (ii-1,jj,dem) ) THEN IF ( dem % mat (ii-1,jj) /= dem % nodata) THEN b = dem % mat (ii-1,jj) ELSE b = dem % mat (ii,jj) END IF ELSE b = dem % mat (ii,jj) END IF !north-east cell IF ( .NOT. IsOutOfGrid (ii-1,jj+1,dem) ) THEN IF ( dem % mat (ii-1,jj+1) /= dem % nodata) THEN c = dem % mat (ii-1,jj+1) ELSE c = dem % mat (ii,jj) END IF ELSE c = dem % mat (ii,jj) END IF !west cell IF ( .NOT. IsOutOfGrid (ii,jj-1,dem) ) THEN IF ( dem % mat (ii,jj-1) /= dem % nodata) THEN d = dem % mat (ii,jj-1) ELSE d = dem % mat (ii,jj) END IF ELSE d = dem % mat (ii,jj) END IF !east cell IF ( .NOT. IsOutOfGrid (ii,jj+1,dem) ) THEN IF ( dem % mat (ii,jj+1) /= dem % nodata) THEN f = dem % mat (ii,jj+1) ELSE f = dem % mat (ii,jj) END IF ELSE f = dem % mat (ii,jj) END IF !south-west cell IF ( .NOT. IsOutOfGrid (ii+1,jj-1,dem) ) THEN IF ( dem % mat (ii+1,jj-1) /= dem % nodata) THEN g = dem % mat (ii+1,jj-1) ELSE g = dem % mat (ii,jj) END IF ELSE g = dem % mat (ii,jj) END IF !south cell IF ( .NOT. IsOutOfGrid (ii+1,jj,dem) ) THEN IF ( dem % mat (ii+1,jj) /= dem % nodata) THEN h = dem % mat (ii+1,jj) ELSE h = dem % mat (ii,jj) END IF ELSE h = dem % mat (ii,jj) END IF !south-east cell IF ( .NOT. IsOutOfGrid (ii+1,jj+1,dem) ) THEN IF ( dem % mat (ii+1,jj+1) /= dem % nodata) THEN i = dem % mat (ii+1,jj+1) ELSE i = dem % mat (ii,jj) END IF ELSE i = dem % mat (ii,jj) END IF !compute the rate of change in the x direction dZdX = ((c + 2.*f + i) - (a + 2.*d + g)) / 8. !compute the rate of change in the y direction dZdY = ((g + 2.*h + i) - (a + 2.*b + c)) / 8. !compute aspect IF (dZdX == 0. .AND. dZdY == 0.) THEN aspect % mat (ii,jj) = -1. CYCLE ELSE aspectIJ = ATAN2 (dZdY, -dZdX) END IF !convert aspect value to compass direction values (0–2Pi radians) IF (aspectIJ < 0) THEN aspect % mat (ii,jj) = pi/2. - aspectIJ ELSE IF (aspectIJ > pi /2.) THEN aspect % mat (ii,jj) = 2.*pi - aspectIJ + pi/2. ELSE aspect % mat (ii,jj) = pi/2. - aspectIJ END IF END IF END DO jdim END DO idim RETURN END SUBROUTINE DeriveAspect !============================================================================== !| Description: ! Calculate curvature according to method implemented in MICROMET ! ! References: ! ! Liston, G.E., Elder, K., A Meteorological Distribution System for ! High-Resolution Terrestrial Modeling, JOURNAL OF HYDROMETEOROLOGY, ! 7, 217-234, 2006. SUBROUTINE CurvatureMicroMet & ! (dem, curve_len_scale, curvature) IMPLICIT NONE !arguments with intent in TYPE(grid_real), INTENT(in) :: dem REAL (KIND = float), INTENT(in) :: curve_len_scale !arguments with intent out TYPE (grid_real), INTENT (out) :: curvature !local declarations: REAL (KIND = float) :: deltaxy REAL (KIND = float) :: z, zw, ze, zs, zn, zsw, zne, znw, zse REAL (KIND = float) :: curve_max INTEGER (KIND = short) :: inc INTEGER (KIND = short) :: i, j !-----------------------------end of declarations------------------------------ !allocate grid CALL NewGrid (curvature, dem) !compute curvature DO i = 1, dem % idim DO j = 1, dem % jdim IF (dem % mat (i,j) /= dem % nodata) THEN !average grid increment deltaxy = CellArea (dem,i,j) ** 0.5 !Convert the length scale to an appropriate grid increment. inc = max(1,nint(curve_len_scale/deltaxy)) !srt z z = dem % mat(i,j) !set zn IF ( .NOT. IsOutOfGrid (i-inc,j,dem) ) THEN IF ( dem % mat (i-inc,j) /= dem % nodata) THEN zn = dem % mat (i-inc,j) ELSE zn = dem % mat (i,j) END IF ELSE zn = dem % mat (i,j) END IF !set zne IF ( .NOT. IsOutOfGrid (i-inc,j+inc,dem) ) THEN IF ( dem % mat (i-inc,j+inc) /= dem % nodata) THEN zne = dem % mat (i-inc,j+inc) ELSE zne = dem % mat (i,j) END IF ELSE zne = dem % mat (i,j) END IF !set ze IF ( .NOT. IsOutOfGrid (i,j+inc,dem) ) THEN IF ( dem % mat (i,j+inc) /= dem % nodata) THEN ze = dem % mat (i,j+inc) ELSE ze = dem % mat (i,j) END IF ELSE ze = dem % mat (i,j) END IF !set zse IF ( .NOT. IsOutOfGrid (i+inc,j+inc,dem) ) THEN IF ( dem % mat (i+inc,j+inc) /= dem % nodata) THEN zse = dem % mat (i+inc,j+inc) ELSE zse = dem % mat (i,j) END IF ELSE zse = dem % mat (i,j) END IF !set zs IF ( .NOT. IsOutOfGrid (i+inc,j,dem) ) THEN IF ( dem % mat (i+inc,j) /= dem % nodata) THEN zs = dem % mat (i+inc,j) ELSE zs = dem % mat (i,j) END IF ELSE zs = dem % mat (i,j) END IF !set zsw IF ( .NOT. IsOutOfGrid (i+inc,j-inc,dem) ) THEN IF ( dem % mat (i+inc,j-inc) /= dem % nodata) THEN zsw = dem % mat (i+inc,j-inc) ELSE zsw = dem % mat (i,j) END IF ELSE zsw = dem % mat (i,j) END IF !set zw IF ( .NOT. IsOutOfGrid (i,j-inc,dem) ) THEN IF ( dem % mat (i,j-inc) /= dem % nodata) THEN zw = dem % mat (i,j-inc) ELSE zw = dem % mat (i,j) END IF ELSE zw = dem % mat (i,j) END IF !set znw IF ( .NOT. IsOutOfGrid (i-inc,j-inc,dem) ) THEN IF ( dem % mat (i-inc,j-inc) /= dem % nodata) THEN znw = dem % mat (i-inc,j-inc) ELSE znw = dem % mat (i,j) END IF ELSE znw = dem % mat (i,j) END IF !curvature curvature % mat(i,j) = (4.0 * z - znw - zse - zsw - zne ) / & (sqrt(2.0) * 16.0 * real(inc) * deltaxy) + & (4.0 * z - zs - zn - ze - zw) / & (16.0 * real(inc) * deltaxy) END IF END DO END DO ! Scale the curvature such that the max abs(curvature) has a value ! of abs(0.5). Include a 1 mm curvature in curve_max to prevent ! divisions by zero in flat terrain where the curvature is zero. curve_max = 0.0 + 0.001 DO i = 1, dem % idim DO j = 1, dem % jdim IF (dem % mat (i,j) /= dem % nodata) THEN curve_max = MAX(curve_max,ABS(curvature%mat(i,j))) END IF END DO END DO DO i = 1, dem % idim DO j = 1, dem % jdim IF (dem % mat (i,j) /= dem % nodata) THEN curvature%mat(i,j) = curvature%mat(i,j) / (2.0 * curve_max) END IF END DO END DO RETURN END SUBROUTINE CurvatureMicroMet !============================================================================== !| Description: ! Calculates the curvature [1/m] of a raster surface, optionally including ! profile and plan curvature. If any neighborhood cells are NoData, they ! are first assigned the value of the center cell. ! It calculates the second derivative value of the input surface on a ! cell-by-cell basis. For each cell, a fourth-order polynomial of the form: !``` ! Z = Ax²y² + Bx²y + Cxy² + Dx² + Ey² + Fxy + Gx + Hy + I !``` ! is fit to a surface composed of a 3x3 window. !``` ! A = [(Z1 + Z3 + Z7 + Z9) / 4 - (Z2 + Z4 + Z6 + Z8) / 2 + Z5] / L4 ! B = [(Z1 + Z3 - Z7 - Z9) /4 - (Z2 - Z8) /2] / L3 ! C = [(-Z1 + Z3 - Z7 + Z9) /4 + (Z4 - Z6)] /2] / L3 ! D = [(Z4 + Z6) /2 - Z5] / L2 ! E = [(Z2 + Z8) /2 - Z5] / L2 ! F = (-Z1 + Z3 + Z7 - Z9) / 4L2 ! G = (-Z4 + Z6) / 2L ! H = (Z2 - Z8) / 2L ! I = Z5 ! ! curvature = -2 (E + D) ! ! profileCurvature = -2 (D G^2 + E H^2 + F G H) / (G^2 + H^2) ! ! planCurvature = -2 (D H^2 + E G^2 - F G H) / (G^2 + H^2) !``` ! Profile curvature is parallel to the direction of the maximum slope. ! A negative value indicates that the surface is upwardly convex at ! that cell. A positive profile indicates that the surface is upwardly ! concave at that cell. A value of zero indicates that the surface ! is linear. Planform curvature (commonly called plan curvature) is ! perpendicular to the direction of the maximum slope. A positive value ! indicates the surface is sidewardly convex at that cell. A negative plan ! indicates the surface is sidewardly concave at that cell. A value of ! zero indicates the surface is linear. Plan curvature relates to the ! convergence and divergence of flow across a surface. ! ! References: ! ! Moore, I. D., R. B. Grayson, and A. R. Landson. 1991. Digital Terrain ! Modelling: A Review of Hydrological, Geomorphological, and Biological ! Applications. Hydrological Processes 5: 3–30. ! ! Zeverbergen, L. W., and C. R. Thorne. 1987. Quantitative Analysis of ! Land Surface Topography. Earth Surface Processes and Landforms 12: 47–56. SUBROUTINE DeriveCurvature & ! (dem, curvature, profileCurvature, planCurvature) IMPLICIT NONE !arguments with intent in TYPE(grid_real), INTENT(in):: dem !arguments with intent out TYPE (grid_real), INTENT (out) :: curvature TYPE (grid_real), OPTIONAL, INTENT (out) :: profileCurvature TYPE (grid_real), OPTIONAL, INTENT (out) :: planCurvature !local variables INTEGER :: i,j REAL (KIND = float) :: D, E, F, G, H REAL (KIND = float) :: z1, z2, z3, z4, z5, z6, z7, z8, z9 REAL (KIND = float) :: L !-----------------------------end of declarations------------------------------ !allocate grids CALL NewGrid (curvature, dem) IF (PRESENT(profileCurvature)) THEN CALL NewGrid (profileCurvature, dem) END IF IF (PRESENT(planCurvature)) THEN CALL NewGrid (planCurvature, dem) END IF !loop through dem grid idim: DO i = 1, dem % idim jdim: DO j = 1, dem % jdim IF (dem % mat(i,j) /= dem % nodata) THEN !set length L = CellArea (dem,i,j) ** 0.5 !set z5 z5 = dem % mat (i,j) !set z2 IF ( .NOT. IsOutOfGrid (i-1,j,dem) ) THEN IF ( dem % mat (i-1,j) /= dem % nodata) THEN z2 = dem % mat (i-1,j) ELSE z2 = dem % mat (i,j) END IF ELSE z2 = dem % mat (i,j) END IF !set z3 IF ( .NOT. IsOutOfGrid (i-1,j+1,dem) ) THEN IF ( dem % mat (i-1,j+1) /= dem % nodata) THEN z3 = dem % mat (i-1,j+1) ELSE z3 = dem % mat (i,j) END IF ELSE z3 = dem % mat (i,j) END IF !set z6 IF ( .NOT. IsOutOfGrid (i,j+1,dem) ) THEN IF ( dem % mat (i,j+1) /= dem % nodata) THEN z6 = dem % mat (i,j+1) ELSE z6 = dem % mat (i,j) END IF ELSE z6 = dem % mat (i,j) END IF !set z9 IF ( .NOT. IsOutOfGrid (i+1,j+1,dem) ) THEN IF ( dem % mat (i+1,j+1) /= dem % nodata) THEN z9 = dem % mat (i+1,j+1) ELSE z9 = dem % mat (i,j) END IF ELSE z9 = dem % mat (i,j) END IF !set z8 IF ( .NOT. IsOutOfGrid (i+1,j,dem) ) THEN IF ( dem % mat (i+1,j) /= dem % nodata) THEN z8 = dem % mat (i+1,j) ELSE z8 = dem % mat (i,j) END IF ELSE z8 = dem % mat (i,j) END IF !set z7 IF ( .NOT. IsOutOfGrid (i+1,j-1,dem) ) THEN IF ( dem % mat (i+1,j-1) /= dem % nodata) THEN z7 = dem % mat (i+1,j-1) ELSE z7 = dem % mat (i,j) END IF ELSE z7 = dem % mat (i,j) END IF !set z4 IF ( .NOT. IsOutOfGrid (i,j-1,dem) ) THEN IF ( dem % mat (i,j-1) /= dem % nodata) THEN z4 = dem % mat (i,j-1) ELSE z4 = dem % mat (i,j) END IF ELSE z4 = dem % mat (i,j) END IF !set z1 IF ( .NOT. IsOutOfGrid (i-1,j-1,dem) ) THEN IF ( dem % mat (i-1,j-1) /= dem % nodata) THEN z1 = dem % mat (i-1,j-1) ELSE z1 = dem % mat (i,j) END IF ELSE z1 = dem % mat (i,j) END IF !compute D, E, F, G, H D = ((z4 + z6) /2. - z5) / L**2 E = ((z2 + z8) /2. - z5) / L**2 F = (-z1 + z3 + z7 - z9) / (4. * L**2) G = (-z4 + z6) / (2. * L) H = (z2 - z8) / (2. * L) curvature % mat (i,j) = -2. * (E + D) IF (PRESENT(profileCurvature)) THEN IF ( (G**2. + H**2.) > 0.) THEN profileCurvature % mat (i,j) = -2. * (D*G**2. + E*H**2. + F*G*H) / (G**2. + H**2.) ELSE profileCurvature % mat (i,j) = 0. END IF END IF IF (PRESENT(planCurvature)) THEN IF ( (G**2. + H**2.) > 0.) THEN planCurvature % mat (i,j) = -2. * (D*H**2. + E*G**2. - F*G*H) / (G**2. + H**2.) ELSE planCurvature % mat (i,j) = 0. END IF END IF END IF END DO jdim END DO idim RETURN END SUBROUTINE DeriveCurvature !============================================================================== !| Description: ! compute map of flow accumulation (m2) ! Input grid: flow direction SUBROUTINE FlowAccumulation & ! (fdir, facc) IMPLICIT NONE !arguments with intent in TYPE(grid_integer), INTENT(in):: fdir !arguments with intent out TYPE (grid_real), INTENT (out) :: facc !local variables INTEGER :: i,j !------------------------------end of declaration ----------------------------- !allocate new flow accumulation grid using flow direction grid as template CALL NewGrid (facc, fdir) !compute grid containing area of each cell DO i = 1, fdir % idim DO j = 1, fdir % jdim IF (fdir % mat (i,j) /= fdir % nodata) THEN facc % mat(i,j) = CellArea (facc, i, j) ELSE facc % mat(i,j) = facc % nodata END IF END DO END DO DO i = 1, fdir % idim DO j = 1, fdir % jdim IF (fdir % mat (i,j) /= fdir % nodata) THEN CALL BasinArea (i, j, fdir, facc % mat(i,j)) ELSE facc % mat(i,j) = facc % nodata END IF END DO END DO END SUBROUTINE FlowAccumulation !============================================================================== !| Description: ! compute basin area (m2) RECURSIVE SUBROUTINE BasinArea & ! (r, c, fdir, area) IMPLICIT NONE TYPE(grid_integer),INTENT(IN) :: fdir INTEGER, INTENt(in) :: r,c REAL (KIND = float), INTENT(inout) :: area !------------------------------end of declaration ----------------------------- IF ( .NOT. IsOutOfGrid(r,c+1,fdir) ) THEN IF(fdir%mat(r,c+1) == W) THEN area = area + CellArea (fdir, r, c+1) CALL BasinArea (r,c+1,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r+1,c+1,fdir) ) THEN IF(fdir%mat(r+1,c+1) == NW ) THEN area = area + CellArea (fdir, r+1, c+1) CALL BasinArea (r+1,c+1,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r+1,c,fdir) ) THEN IF(fdir%mat(r+1,c) == N ) THEN area = area + CellArea (fdir, r+1, c) CALL BasinArea (r+1,c,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r+1,c-1,fdir) ) THEN IF(fdir%mat(r+1,c-1) == NE ) THEN area = area + CellArea (fdir, r+1, c-1) CALL BasinArea (r+1,c-1,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r,c-1,fdir) ) THEN IF(fdir%mat(r,c-1) == E) THEN area = area + CellArea (fdir, r, c-1) CALL BasinArea (r,c-1,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r-1,c-1,fdir) ) THEN IF(fdir%mat(r-1,c-1) == SE ) THEN area = area + CellArea (fdir, r-1, c-1) CALL BasinArea (r-1,c-1,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r-1,c,fdir) ) THEN IF(fdir%mat(r-1,c) == S ) THEN area = area + CellArea (fdir, r-1, c) CALL BasinArea (r-1,c,fdir,area) END IF END IF IF ( .NOT. IsOutOfGrid(r-1,c+1,fdir) ) THEN IF(fdir%mat(r-1,c+1) == SW ) THEN area = area + CellArea (fdir, r-1, c+1) CALL BasinArea (r-1,c+1,fdir,area) END IF END IF END SUBROUTINE BasinArea !============================================================================== !| Description: ! find the cells that are springs, defined as those cells that have not ! any other cells upstream FUNCTION CellIsSpring & ! (row, col, flowDir) & ! RESULT (spring) IMPLICIT NONE ! Arguments with intent(in): INTEGER, INTENT(in) :: row, col TYPE (grid_integer), INTENT(in) :: flowDir ! Local variables: LOGICAL :: spring !------------end of declaration------------------------------------------------ spring = .TRUE. IF (flowDir % mat(row,col) == flowDir % nodata) THEN spring = .FALSE. RETURN END IF IF(.NOT. IsOutOfGrid(row,col+1,flowDir) ) THEN IF(flowDir%mat(row,col+1) == W ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row+1,col+1,flowDir) ) THEN IF(flowDir%mat(row+1,col+1) == NW ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row+1,col,flowDir) ) THEN IF(flowDir%mat(row+1,col) == N ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row+1,col-1,flowDir) ) THEN IF(flowDir%mat(row+1,col-1) == NE ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row,col-1,flowDir) ) THEN IF(flowDir%mat(row,col-1) == E ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row-1,col-1,flowDir) ) THEN IF(flowDir%mat(row-1,col-1) == SE ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row-1,col,flowDir) ) THEN IF(flowDir%mat(row-1,col) == S ) THEN spring = .FALSE. RETURN ENDIF ENDIF IF(.NOT. IsOutOfGrid(row-1,col+1,flowDir) ) THEN IF(flowDir%mat(row-1,col+1) == SW ) THEN spring = .FALSE. RETURN ENDIF ENDIF END FUNCTION CellIsSpring !============================================================================== !| Description: ! compute mask of river basin given map of flow direction and the coordinate ! of the outlet point SUBROUTINE BasinDelineate & ! (fdir,x,y, mask) IMPLICIT NONE !arguments with intent in TYPE(grid_integer), INTENT(in):: fdir REAL (KIND = float), INTENT(in) :: x, y !!coordinate of outlet !arguments with intent inout: TYPE(grid_integer), INTENT(inout) :: mask !local variables INTEGER :: i,j LOGICAL :: check !------------------------------end of declaration ----------------------------- IF (.NOT. ALLOCATED (mask % mat)) THEN !allocate new grid CALL NewGrid (mask, fdir) END IF mask % mat = mask % nodata !compute row and column of outlet CALL GetIJ (x, y, fdir, i, j, check) !compute basin mask CALL BasinMask(mask,fdir,i,j) END SUBROUTINE BasinDelineate !============================================================================== !| Description: ! search for cells included in the river basin RECURSIVE SUBROUTINE BasinMask & ! (basin, fdir, r, c) IMPLICIT NONE TYPE(grid_integer),INTENT(IN) :: fdir TYPE(grid_integer),INTENT(INOUT) :: basin INTEGER, INTENT(in) :: r,c !------------------------------end of declaration ----------------------------- IF ( .NOT. IsOutOfGrid(r,c+1,fdir) ) THEN IF(fdir%mat(r,c+1)==W.AND. basin%mat(r,c+1)/=1) THEN basin%mat(r,c+1) = 1 CALL BasinMask(basin,fdir,r,c+1) END IF END IF IF ( .NOT. IsOutOfGrid(r+1,c+1,fdir) ) THEN IF(fdir%mat(r+1,c+1)==NW .AND. basin%mat(r+1,c+1)/=1) THEN basin%mat(r+1,c+1) = 1 CALL BasinMask(basin,fdir,r+1,c+1) END IF END IF IF ( .NOT. IsOutOfGrid(r+1,c,fdir) ) THEN IF(fdir%mat(r+1,c)==N .AND. basin%mat(r+1,c)/=1) THEN basin%mat(r+1,c) = 1 CALL BasinMask(basin,fdir,r+1,c) END IF END IF IF ( .NOT. IsOutOfGrid(r+1,c-1,fdir) ) THEN IF(fdir%mat(r+1,c-1)==NE .AND. basin%mat(r+1,c-1)/=1) THEN basin%mat(r+1,c-1) = 1 CALL BasinMask(basin,fdir,r+1,c-1) END IF END IF IF ( .NOT. IsOutOfGrid(r,c-1,fdir) ) THEN IF(fdir%mat(r,c-1)==E .AND. basin%mat(r,c-1)/=1) THEN basin%mat(r,c-1) = 1 CALL BasinMask(basin,fdir,r,c-1) END IF END IF IF ( .NOT. IsOutOfGrid(r-1,c-1,fdir) ) THEN IF(fdir%mat(r-1,c-1)==SE .AND. basin%mat(r-1,c-1)/=1) THEN basin%mat(r-1,c-1) = 1 CALL BasinMask(basin,fdir,r-1,c-1) END IF END IF IF ( .NOT. IsOutOfGrid(r-1,c,fdir) ) THEN IF(fdir%mat(r-1,c)==S .AND. basin%mat(r-1,c)/=1) THEN basin%mat(r-1,c) = 1 CALL BasinMask(basin,fdir,r-1,c) END IF END IF IF ( .NOT. IsOutOfGrid(r-1,c+1,fdir) ) THEN IF(fdir%mat(r-1,c+1)==SW .AND. basin%mat(r-1,c+1)/=1) THEN basin%mat(r-1,c+1) = 1 CALL BasinMask(basin,fdir,r-1,c+1) END IF END IF END SUBROUTINE BasinMask !============================================================================== !| Description: ! compute longest flow length (m) FUNCTION LongestFlowLength & ! (fdir, x, y) & ! RESULT (lmax) IMPLICIT NONE !Arguments with intent (in) TYPE(grid_integer),INTENT(IN) :: fdir REAL (KIND = float), INTENT(in) :: x, y !local declarations REAL (KIND = float) :: lmax REAL (KIND = float) :: length REAL (KIND = float) :: xu, yu, xd, yd TYPE(grid_integer) :: basin INTEGER (KIND = short) :: row, col !!current cell INTEGER (KIND = short) :: iDown, jDown !!downstream cell INTEGER (KIND = short) :: i, j LOGICAL :: outlet !------------------------------end of declaration ----------------------------- !delineate river basin CALL BasinDelineate (fdir, x, y, basin) !overlay flow direction map CALL GridResample (fdir, basin) point1 % system = basin % grid_mapping point2 % system = basin % grid_mapping lmax = 0. !loop trough basin DO j = 1,basin % jdim DO i = 1,basin % idim IF (basin % mat (i,j) /= basin % nodata) THEN IF(CellIsSpring(i,j,basin)) THEN !found a spring length = 0. !follow the reach till basin outlet row = i col = j outlet = .FALSE. DO WHILE (.NOT. outlet) ! follow the reach till the basin outlet CALL DownstreamCell(row, col, basin%mat(row,col), iDown, jDown) CALL GetXY (row,col,basin,xu,yu) CALL GetXY (iDown,jDown,basin,xd,yd) point1 % northing = yu point1 % easting = xu point2 % northing = yd point2 % easting = xd length = length + Distance(point1,point2) outlet = CheckOutlet(row,col,iDown,jDown,basin) IF (outlet) THEN IF (length > lmax) THEN lmax = length END IF END IF !loop row = iDown col = jDown END DO ENDIF END IF ENDDO ENDDO RETURN END FUNCTION LongestFlowLength !============================================================================== !| Description: ! compute slope of longest flow path (m/m) FUNCTION LongestFlowPathSlope & ! (fdir, dem, x, y, outsection, headsection) & ! RESULT (slope) IMPLICIT NONE !Arguments with intent (in) TYPE(grid_integer),INTENT(IN) :: fdir !flow direction TYPE(grid_real),INTENT(IN) :: dem !digital elevation model REAL (KIND = float), INTENT(in) :: x, y !Arguments with intent out TYPE (Coordinate), OPTIONAL, INTENT (OUT) :: outsection !!coordinate of outlet section TYPE (Coordinate), OPTIONAL, INTENT (OUT) :: headsection !!coordinate of head section !local declarations REAL (KIND = float) :: slope REAL (KIND = float) :: lmax REAL (KIND = float) :: length REAL (KIND = float) :: xu, yu, xd, yd TYPE(grid_integer) :: basin INTEGER (KIND = short) :: row, col !!current cell INTEGER (KIND = short) :: iDown, jDown !!downstream cell INTEGER (KIND = short) :: iSpring, jSpring !!coordinate of the spring INTEGER (KIND = short) :: iHead, jHead !!coordinate of the head of the longest flow path INTEGER (KIND = short) :: iOutlet, jOutlet !!coordinate of the basin outlet INTEGER (KIND = short) :: i, j LOGICAL :: outlet !------------------------------end of declaration ----------------------------- !delineate river basin CALL BasinDelineate (fdir, x, y, basin) !overlay flow direction map CALL GridResample (fdir, basin) !find local coordinate of basin outlet CALL GetIJ (x, y, fdir, iOutlet, jOutlet) point1 % system = basin % grid_mapping point2 % system = basin % grid_mapping lmax = 0. iHead = 0 jHead = 0 !loop trough basin DO j = 1,basin % jdim DO i = 1,basin % idim IF (basin % mat (i,j) /= basin % nodata) THEN IF(CellIsSpring(i,j,basin)) THEN !found a spring length = 0. iSpring = i jSpring = j !follow the reach till basin outlet row = i col = j outlet = .FALSE. DO WHILE (.NOT. outlet) ! follow the reach till the basin outlet CALL DownstreamCell(row, col, basin%mat(row,col), iDown, jDown) CALL GetXY (row,col,basin,xu,yu) CALL GetXY (iDown,jDown,basin,xd,yd) point1 % northing = yu point1 % easting = xu point2 % northing = yd point2 % easting = xd length = length + Distance(point1,point2) outlet = CheckOutlet(row,col,iDown,jDown,basin) IF (outlet) THEN IF (length > lmax) THEN !update lmax lmax = length !update head iHead = iSpring jHead = jSpring END IF END IF !loop row = iDown col = jDown END DO ENDIF END IF ENDDO ENDDO !compute slope (m/m) slope = ( dem % mat(iHead,jHead) - dem % mat (iOutlet,jOutlet) ) / lmax IF (PRESENT(outsection)) THEN outsection % system = basin % grid_mapping CALL GetXY (iOutlet, jOutlet, basin, outsection % easting, outsection % northing) END IF IF (PRESENT(headsection)) THEN headsection % system = basin % grid_mapping CALL GetXY (iHead, jHead, basin, headsection % easting, headsection % northing) END IF RETURN END FUNCTION LongestFlowPathSlope !============================================================================== !| Description: ! Drainage Density (1/m) of a basin is the total line length of all streams ! divided by basin area. If flow accumulation is not given, it is computed ! from flow direction. If coordinate of closing section is not given ! drainage density is computed for the entire grid. FUNCTION DrainageDensity & ! (channel, fdir, x, y, flowAcc) & ! RESULT (dd) IMPLICIT NONE !Arguments with intent (in) TYPE(grid_integer),INTENT(IN) :: channel !!mask of channel cells. NoData for non channel cells TYPE(grid_integer),INTENT(IN) :: fdir !!flow direction REAL (KIND = float), INTENT(IN),OPTIONAL :: x, y !!coordinate of basin closing section TYPE (grid_real),INTENT(IN),OPTIONAL :: flowAcc !!flow accumulation !local declarations REAL (KIND = float) :: dd TYPE (grid_real) :: facc TYPE (grid_integer) :: mask INTEGER (KIND = short) :: i,j INTEGER (KIND = short) :: is, js !!coordinate of downstream cell in local reference system REAL (KIND = float) :: length !!total river length REAL (KIND = float) :: cLength !!cell length LOGICAL :: check !------------------------------end of declaration ----------------------------- !set flowaccumulation IF (PRESENT (flowAcc)) THEN CALL NewGrid (facc, fdir) facc = flowAcc ELSE CALL FlowAccumulation (fdir, facc) END IF IF (PRESENT(x) .AND. PRESENT(y)) THEN !delineate watershed mask CALL BasinDelineate (fdir,x,y, mask) ELSE !drainage density is computed on entire grid CALL NewGrid (mask, fdir) END IF !compute total channel length [m] length = 0. DO i = 1, mask % idim DO j = 1, mask % jdim IF (mask % mat (i,j) /= mask % nodata) THEN IF (channel % mat (i,j) /= channel % nodata) THEN CALL DownstreamCell (i, j, fdir % mat(i,j), is, js, dx = cLength, grid = fdir) length = length + cLength END IF END IF END DO END DO !compute drainage density [1/m] CALL GetIJ (x, y, fdir, i, j, check) IF (.NOT. check) THEN CALL Catch ('error', 'Morphology', 'Drainage density computation: ', & argument = 'closing section outside grid dimension') END IF IF (channel % mat(i,j) == channel % nodata) THEN CALL Catch ('error', 'Morphology', 'Drainage density computation: ', & argument = 'closing section on hillslope') END IF dd = length / facc % mat(i,j) !deallocate CALL GridDestroy (facc) CALL GridDestroy (mask) RETURN END FUNCTION DrainageDensity !============================================================================== !| Description: ! Topographic index `TI` takes into account both a local slope geometry ! and site location in the landscape combining data on slope steepness `S` ! and specific catchment area `As`: !``` ! TI = ln (As/S) !``` ! For grid structures the specific catchment area can be obtained by ! dividing the contributing area of the cell by the effective contour length. ! This is the length of contour line within the grid cell over which flow ! can pass. It is calculated as the length of the line through the grid ! cell centre and perpendicular to the aspect direction. ! ! References: ! ! P. J. J. DESMET & G. GOVERS (1996): Comparison of routing algorithms ! for digital elevation models and their implications for predicting ! ephemeral gullies, International Journal of Geographical ! Information Systems, 10:3, 311-331 SUBROUTINE TopographicIndex & ! (dem, fdir, facc, xcoord, ycoord, tivalue, tigrid) IMPLICIT NONE !Argumentd with intent(in): TYPE (grid_real), INTENT(IN) :: dem !!digital elevation model TYPE (grid_integer), INTENT(IN) :: fdir !!flow direction TYPE (grid_real), OPTIONAL, INTENT(IN) :: facc !!flow accumulation in cell unit REAL (KIND = float) , OPTIONAL, INTENT(IN) :: xcoord !! x coordinate of cell for computing index REAL (KIND = float) , OPTIONAL, INTENT(IN) :: ycoord !! y coordinate of cell for computing index !Arguments with intent(out): REAL (KIND = float) , OPTIONAL, INTENT(OUT) :: tivalue !! index of cell with coordinates (x,y) TYPE (grid_real), OPTIONAL, INTENT(OUT) :: tigrid !!topographic index !local declarations TYPE (grid_real) :: area !! contributing area (m^2) TYPE (grid_real) :: aspect !!aspect (rad) TYPE (grid_real) :: slope !!slope (rad) REAL (KIND = float) :: as !! specific catchment area (m^2 / m) REAL (KIND = float) :: x !!factor converting cellsize along direction... !!...perpendicular to the aspect direction REAL (KIND = float) :: alpha !! local aspect (rad) REAL (KIND = float) :: cellsize !!cell size (m) INTEGER (KIND = short) :: i, j !------------------------------end of declaration ----------------------------- !allocate topographic index grid CALL NewGrid (tigrid, dem) IF ( .NOT. PRESENT(facc) ) THEN !compute contributing area CALL FlowAccumulation (fdir, area) ELSE CALL NewGrid (area, fdir) DO i = 1, area % idim DO j = 1, area % jdim IF (area % mat (i,j) /= area % nodata ) THEN area % mat (i,j) = facc % mat (i,j) ! (* CellArea (facc, i, j)) END IF END DO END DO END IF !compute aspect CALL DeriveAspect (dem, aspect) !compute slope CALL DeriveSlope (dem, slope) !if present tivalue compute topographic index only for cell with coordinate xcoord, ycoord IF (PRESENT (tivalue) ) THEN !check wether coordinates of cell is given IF ( PRESENT (xcoord) .AND. PRESENT (ycoord) ) THEN CALL GetIJ (xcoord, ycoord, fdir, i, j) alpha = aspect % mat (i,j) x = 1. / MAX(ABS(SIN(alpha)),ABS(COS(alpha))) cellsize = ( CellArea (fdir, i, j) ) ** 0.5 ![m] as = area % mat (i,j) / (cellsize * x) IF ( TAN (slope % mat (i,j)) == 0.) THEN tivalue = LOG (as / 0.0001 ) ELSE tivalue = LOG (as / TAN (slope % mat (i,j))) END IF ELSE CALL Catch ('error', 'Morphology', 'missing coordinate for compuing topographic index of given cell') END IF END IF !If present tigrid, compute topographic index for every cell of the grid IF (PRESENT (tigrid) ) THEN DO i = 1, tigrid % idim DO j = 1, tigrid % jdim IF (tigrid % mat (i,j) /= tigrid % nodata) THEN alpha = aspect % mat (i,j) x = 1. / MAX(ABS(SIN(alpha)),ABS(COS(alpha))) cellsize = ( CellArea (tigrid, i, j) ) ** 0.5 ![m] as = area % mat (i,j) / (cellsize * x) IF ( TAN (slope % mat (i,j)) == 0.) THEN tigrid % mat (i,j) = LOG (as / 0.0001) ELSE tigrid % mat (i,j) = LOG (as / TAN (slope % mat (i,j))) END IF END IF END DO END DO END IF !purge CALL GridDestroy (area) CALL GridDestroy (aspect) CALL GridDestroy (slope) RETURN END SUBROUTINE TopographicIndex !============================================================================== !| Description: ! find the cells that are springs, defined as those cells that have not ! any other cells upstream FUNCTION Perimeter & ! (grid) & ! RESULT (p) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: grid ! Local variables: REAL (KIND = float) :: p TYPE (grid_integer) :: border INTEGER :: i,j !------------end of declaration------------------------------------------------ CALL ExtractBorder (grid, border, cardinal = .TRUE.) p = 0. DO i = 1, border % idim DO j = 1, border % jdim IF (border % mat (i,j) /= border % nodata ) THEN p = p + CellLength (border, i,j) END IF END DO END DO !increase a default value of 9% to consider diagonal length p = p * 1.09 CALL GridDestroy (border) RETURN END FUNCTION Perimeter END MODULE Morphology